/*
 * This file contains all functions in the course Comp. Phys., including:
 ************************************************************************
 *   * job *			* algorithm *		* function *
 ************************************************************************
 * derivative of a 1-dim function:		 	Richardson_deriv
 * interpolation/extrapolation:				Nevill
 * numerical integration:	midpoint algorithm	integ_midpoint
 * 				trapzoid		integ_trapzoid
 * 				Simpson			Simpson
 * 				SelfAdaptive Simpson	SelfAdapt_Simpson
 * 				Gauss-Legendre		GauLeg
 * 							GauLegInt
 * nonlinear equation:		bisection algorithm	bisection
 * 				Newton method		Newton_findroot
 * 				secant method		SecantMethod
 * 1-dim optimization					findmin
 * equation set			Gaussian elimination	GaussLinear
 * 				LU decomposition	ludcmp
 * 							ludcmp_solve
 * 				tridiagonal case	TridiagLinear (not rebust yet)
 * diagonalization		jacobi algorithm	jacobi
 */

#define max_iteration 1E6

/*
 * Richardson算法：
 * 首先计算
 *      D^k_0 = (f(x+h/2^k) - f(x-h/2^k))/ ( h/2^{k-1} ), k=0,1,...,n
 * 然后根据递推公式
 *      D^k_m = D^{k+1}_{m-1} + \frac{ D^{k+1}_{m-1} - D^k_{m-1} }{ 4^m -1 }
 * 逐步推得　D^0_n 作为返回值，其误差为　O(h^{2(n+1)}
 *
 * 例：
 *
 * D^0_0
 * D^1_0        D^0_1
 * D^2_0        D^1_1   D^0_2
 *
 * 我们定义一个　3x3　的矩阵 deriv，然后逐列计算矩阵元，最后返回　D^0_2, 即 deriv[2][2] 的值。
 *
 * n 不能是负数，否则程序终止
 */
double Richardson_deriv(double x, double h, int n, double (*func)(double) ){

        double coef, y;

        double ** deriv = new double * [ n+1 ];//开辟动态内存，得到新的 (n+1)*(n+1) 矩阵
        for(int i=0;i<n+1;i++) deriv[i] = new double [n+1];

        for(int i=0;i<=n;i++){ // 赋值第一列 D^k_0, k=0,1,...,n
                coef = pow(0.5,i);// 1/(2^i)
                deriv[i][0] = ( func(x + h*coef ) - func(x - h*coef) ) / (2* h* coef);
        }
        for(int i=1;i<=n;i++){ // 赋值第　1,...,n 列
                for(int j=i;j<=n;j++){ // j 是行数，i是列数
                        deriv[j][i] = deriv[j][i-1] + (deriv[j][i-1] - deriv[j-1][i-1])/ (pow(4.0,i) -1 );//　这一行与递推公式对应
                }
        }
        //cout<<" error is O(h^"<<2*n+2<<")\n";
        y = deriv[n][n];
        for(int i=0;i<=n;i++)delete [] deriv[i];//释放动态内存
        delete [] deriv;

        return y;
}

/*
 * Nevill 算法
 *
 * 给定 n 个数据点　(x_data[], y_data[])，指定内插/外推点位置　x，Nevill()函数自动进行多项式拟合，得>到　x　处插值结果 y，并计算误差 dy
 *
 * 如果数据点的　x　坐标太靠近，误差可能会比较大，所以程序中断。
 *
 */
double Nevill(int n, double * x_data, double * y_data, double x, double &y, double & dy){

        double *c = new double [n];
        double *d = new double [n];

        int i, m, ns=0;
        double dis=fabs(x-x_data[0]), dif, den;
        for(int i=0;i<n;i++){
                if( dif=fabs(x-x_data[i]) < dis ){
			// 找到离x最近的数据点x[ns]，从它的函数值开始逐步推测y
                        dis = dif;
                        ns = i;
                }
                c[i] = y_data[i];//赋初值　c,d,　可以证明，这个初值带入公式，确然得到第二列　P_{01}, P_{12}, ...
                d[i] = y_data[i];
        }
        y = y_data[ns];//y的初始值

        for(m=1;m<n;m++){
                for(i=0;i<n-m;i++){
                        den = x_data[i] - x_data[i+m];
                        if( fabs(den) < 1E-9 ){
                                cout<<" two data points' x coordinates are too close.\n";
                                exit(1);
                        }
                        den = ( c[i+1] - d[i] )/ den;

                        d[i] = (x_data[i+m] - x) * den;//迭代
                        c[i] = (x_data[i] - x) * den;
                }
                y += ( dy = (2*ns < n-m ? c[ns]: d[--ns]) );// 这一行保证推测路径最接近直线。 ns < (n-m)/2，说明 ns 在上半区，由　c路径 （斜向下），否则由　d路径（斜向上）
        }
        delete [] c;
        delete [] d;
}
/*
 * integ_midpoint: 中值法积分
 * 给定区间[a,b]，分割为n份，被积函数为func，返回积分值
 */
double integ_midpoint(double a, double b, int n, double (*func)(double)){

        int i;
        double y=0;
        double h=(b-a)/n;
        for(i=0;i<n;i++){
                y += func( a + (i+0.5)*h );
        }
        y *= h;
        return y;
}
/*
 * integ_trapzoid: 梯形法积分
 * 给定区间[a,b]，分割为n份，被积函数为func，返回积分值
 */
double integ_trapzoid(double a, double b, int n, double (*func)(double)){

        int i;
        double y=func(a)/2;
        double h=(b-a)/n;
        for(i=1;i<=n-1;i++){
                y += func( a + i*h );
        }
        y += func(b)/2;
        y *= h;
        return y;
}

/*
 * Simpson：辛普生积分
 * 给定积分区间[a,b]，区间分割为n份，被积函数为func，返回积分值
 */
double Simpson(double a, double b, int n, double (*func)(double)){

        double y1,y2,y,h;
        int i;

        h = (b-a)/n;//步长

        y1 = 0;
        for(i=0;i<n;i++){
                y1 += func(a+ (i+0.5)*h);
        }
        y2 = 0;
        for(i=1;i<n;i++){
                y2 += func(a+ i*h);
        }
        y = 4*y1 + 2*y2 + func(a) + func(b);
        y *= h/6;
        return y;
}

/*
 * SelfAdapt_Simpson uses self adapted algorithm to calculate integration of func() over [a,b].
 * It automatically breaks [a,b] into unequal pieces, according to the local situaion of func().
 * the error is controled to be less than epsilon
 */
double SelfAdapt_Simpson(double a, double b, double (*func)(double), double epsilon){

        double y1 = Simpson(a, (a+b)/2, 1, func);
        double y2 = Simpson((a+b)/2, b, 1, func);
        double y3 = Simpson(a, b, 1, func);
        if( fabs(y1+y2-y3) < 15*epsilon ){
                return y3;
        }
        else{
                return SelfAdapt_Simpson(a, (a+b)/2, func, epsilon/2)
                        +
                        SelfAdapt_Simpson((a+b)/2, b, func, epsilon/2);
        }
}

/*
 * Given [a,b], returns abscissas x[] and weights w[], based on zeros of P_n(x)
 * Modified from Numerical Recipes
 */
void GauLeg(double a, double b, int n, double *x, double *w){

        double eps=1E-14;
        int m,j,i;
        double z1,z,xm,xl,pp,p3,p2,p1;

        m=(n+1)/2;
        xm=0.5*(a+b);
        xl=0.5*(b-a);
        for(i=0;i<m;i++){
                z=cos(M_PI*(i+0.75)/(n+0.5));//starting guess
                do{                     //Newton's method to get the ith zero
                        p1=1.0;         // P_0
                        p2=0.0;
                        for(j=0;j<n;j++){
                                p3=p2;  // P_{j-1}
                                p2=p1;  // P_{j}
                                p1= z*p2 + (z*p2-p3)*j/(j+1);
                                // P_{j+1},based on recurrence relation:
                                // (j+1)P_{j+1} = (2j+1)P_j - j P_{j-1}
                        }
                        pp=n*(z*p1-p2)/(z*z-1); // P'n(x_i), based on
                                    // (1-x^2)P'_n(x) = n P_{n-1}(x) - nxP_n (x)
                        z1=z;
                        z=z1-p1/pp; // Newton's method: s = x - f(x)/f'(x)
                }while(fabs(z-z1)>eps);
                x[i]=xm-xl*z;
                x[n-1-i]=xm+xl*z;
                w[i] = 2.0*xl/((1.0-z*z)*pp*pp); // formula of weights:
                                    // w[i] = 2/(1-x_j^2)/P'n(x_j)^2 for [-1,1]
                                    // multiplied by xl, for [a,b]
                w[n-1-i]=w[i];
        }
}

/*
 * Gauss-Legendre Integral
 * Order-n Gauss-Legendre integral, returns integral of func on [a,b]
 */
double GauLegInt(double a, double b, double (*func)(double), int n){

        if(n<1){
                cout<<"error: n<1 in GauLegInt(a,b,func,n).\n";
                exit(1);
        }
        int i;
        double y;

        double *x = new double [n];
        double *w = new double [n];

        GauLeg(a,b,n,x,w);//get the abscissas and weights

        y=0;
        for(i=0;i<n;i++){
                y += w[i] * func(x[i]);
        }

        delete [] x;
        delete [] w;

        return y;
}

/*
 * findroot: 在[a,b]区间内寻找 func(x)=0　的根，并返回根的值
 * 算法：二分法
 */
double bisection(double (*func)(double), double a, double b, double precision){

        double ya=func(a), yb=func(b);
        double c=(a+b)/2, yc;
        int i;

        if( ya* yb > 0){//如果 func(a) * func(b) >0, 则[a,b]区间内不一定有解，拒绝任务
                cout<<"error: invalid input [a,b] for findroot: f(a)*f(b)>0 "<<endl;
                exit(1);
        }

        for(i=0;i<max_iteration;i++){

                c = (a+b)/2;//区间中点
                yc = func(c);

                if( yc * ya <= 0 ){//如果 yc*ya<=0，说明[a,c]中有一个根
                        b=c;
                }
                else{//如果 yc*ya>0，说明[c,b]中有一个根
                        a=c;
                }
                if( fabs(a-b) < precision ){
			cout<<"\n iteration rounds = "<<i+1<<endl;
                        return (a+b)/2;
                }
        }
        cout<<"bisection method: failed to find a root, after "<<max_iteration<<" iterations.\n";
        exit(1);
}

/*
 * Newton's method: given a function (*func), which gives value and also the 1st derivative,
 * a suspicious zone [a,b], it finds a root in the zone within precision, and returns the root.
 */
double Newton_findroot(void (*func)(double, double &, double &), double a, double b, double precision){

        double x, dx, y, dy;
        int i;

        x = (a+b)/2;

        for(i=0;i<max_iteration;i++){
                func(x, y, dy);
                if( fabs(dy)<1E-9 ){
                        cout<<" error: f'=0 in Newton's root-finding method. \n";
                        exit(1);
                }
                dx = y/dy;
                x -= dx;
                if( (x-a)*(x-b) >0 ){
                        cout<<" error: x gets out of the given zone [a,b] in Newton's method.\n";
                        exit(1);
                }
                if( fabs(dx) < precision ){
			cout<<"\n iteration rounds = "<<i+1<<endl;
                        return x;
                }
        }
        cout<<" Newton's method: after "<<max_iteration<<" steps, failed to find a root.\n";
        exit(1);
}


/*
 * SecantMethod: given function func, suspicious zone [a,b], returns a root of func(x)==0
 */
double SecantMethod(double (*func)(double), double a, double b, double precision ){

        double c, fc, d, fd, e, fe;
        int i;

        c = (a+b)/2;//set starting points
        fc = func(c);
        d = (c+b)/2;
        fd = func(d);

        for(i=0;i<max_iteration;i++){
                e = d - (d-c)/(fd-fc)*fd;
                if( (e-a)*(e-b) >0 ){
                        cout<<"SecantMethod: error! iteration gets out of the given zone [a,b].\n";
                        exit(1);
                }
                fe = func(e);
                if( fabs(e-d)<precision && fabs(fe)<precision ){
			cout<<"\n iteration rounds = "<<i+1<<endl;
                        return e;
                }
                c=d;
                fc=fd;
                d=e;
                fd=fe;
        }

        cout<<"After "<<max_iteration<<" iterations, failed to find one root"<<endl;
        exit(1);
}

/*
 * StewensonMethod: given function func, suspicious zone [a,b], returns a root of func(x)==0
 */
double StewensonMethod(double (*func)(double), double a, double b, double precision ){

        double d, fd, e, fe;
        int i;

        d = (a+b)/2;
        fd = func(d);

        for(i=0;i<max_iteration;i++){
		if( fabs(fd)< 1E-15 ){
			cout<<"\n iteration rounds = "<<i+1<<endl;
			return d;
		}
                e = d - fd* fd/(fd - func(d - fd));
                if( (e-a)*(e-b) >0 ){
                        cout<<"StewensonMethod: error! iteration gets out of the given zone [a,b].\n";
                        exit(1);
                }
                fe = func(e);
                if( fabs(e-d)<precision ){
			cout<<"\n iteration rounds = "<<i+1<<endl;
                        return e;
                }
                d=e;
                fd=fe;
        }

        cout<<"After "<<max_iteration<<" iterations, failed to find one root"<<endl;
        exit(1);
}


/*
 * findmin: 给定函数func, 初始区间[a,b]，在精度要求precision下，寻找一个极小值点，并返回该点坐标
 */
double findmin(double (*func)(double), double a, double b, double precision){

        double temp,c,d;
        int i;

        if( func(a)<func(b) ){//寻找极小值的方向为：a->b->
                temp=a;
                a=b;
                b=temp;
        }

        for(i=0;i<max_iteration;i++){//寻找嫌疑区间
                c= b + 1.618*(b-a);
                if( func(c) > func(b) )break;
                a=b;
                b=c;
        }
        if(i==max_iteration){//10000次迭代以后，没有找到嫌疑区间，报错
                cout<<"after "<<max_iteration<<" iterations, didn't find a zone for fmin"<<endl;
                exit(1);
        }

        for(i=0;i<max_iteration;i++){
                if( fabs(a-c)<precision )break;
                d=0.38197*a + 0.61803*c;
                if( (d-a)*(d-b)<0){//如果d在a,b之间
                        if( func(d) < func(b) ){
                                c=b;
                                b=d;
                        }
                        else{
                                a=d;
                        }
                }
                else if( (d-a)*(d-b)>0 ){//如果d在b,c区间
                        if( func(d) < func(b) ){
                                a=b;
                                b=d;
                        }
                        else{
                                c=d;
                        }
                }
                //下面四行帮忙查错：观察每次迭代的结果
                //cout<<"-----------------------------------"<<endl;
                //cout<<"a="<<a<<", f(a)="<<func(a)<<endl;
                //cout<<"b="<<a<<", f(b)="<<func(b)<<endl;
                //cout<<"c="<<a<<", f(c)="<<func(c)<<endl;
        }
        return (a+c)/2;
}


/*
 * GaussLinear: solves n-dim linear equations
 * The solution is stored in x
 * Permutation of rows is activated, when the principle element is smaller than min
 * If A is singular, a permutation will fail and function exit(1), signals error
 */
void GaussLinear(int n, double **A, double *x, double min){

        int i,j,k;
        double y;
        for(i=0;i<n;i++){
                if( fabs( A[i][i] ) < min ){// permutation is activated
                        for(j=i+1;j<n;j++){
                                if( fabs( A[j][i] ) > min ){
                                        for(k=i;k<n;k++){
                                                y = A[i][k];
                                                A[i][k] = A[j][k];
                                                A[j][k] = y;
                                        }
                                        break;
                                }
                        }
                        if(j==n){// didn't find any element larger than min
                                cerr<<" A is singular in GaussLinear(n,A,x,min).\n";
                                exit(1);
                        }
                }
                for(j=i+1;j<n;j++){// row operations
                        y = - A[j][i] / A[i][i];
                        for(k=0;k<n+1;k++){
                                A[j][k] = y * A[i][k] + A[j][k];
                        }
                }
                // for test only
                /*
                cout<<"i="<<i<<"\t A:"<<endl;
                for(j=0;j<n;j++){
                        for(k=0;k<n+1;k++){
                                cout<<A[j][k]<<",";
                        }
                        cout<<endl;
                }
                */
        }
        for(j=n-1;j>=0;j--){
                y=0;
                for(k=j+1;k<n;k++) y += A[j][k]*x[k];
                x[j] = 1/A[j][j] * ( A[j][n] - y );
        }
}

/*
 * ludcmp(): 	LU decomposition of a non-singular matrix A
 * The L matrix has L_{ii} = 1, and other elements are written into lower-left part of a, while the U matrix is written into the upper-right part of a.
 * Note that ludcmp does not rule out all singular cases. Some singular A can also be successfully decomposed by ludcmp.
 * n:		dimension
 * a:		n x n matrix
 * pivot:	records pivotting
 * sign:	1/-1 for even/odd times of pivotting
 */
void ludcmp(int n, double **a, int *pivot, int & sign){

	int i,j,k,imax;
	double big, y;
	double *max = new double [n];
	for(i=0;i<n;i++){
		big = 0;
		for(j=0;j<n;j++){
			if( fabs(a[i][j]) > big ){
				big = fabs( a[i][j] );
			}
		}
		if( big < 1E-9 ){
			cout<<"a is singular in ludcmp(n,a,pivot,sign).\n";
			exit(1);
		}
		max[i] = big;
	}
	//cout<<"max:"; for(i=0;i<n;i++)cout<<max[i]<<","; cout<<endl;

	sign =1;
	for(j=0;j<n;j++){
		for(i=0;i<=j;i++){
			y=0;
			for(k=0;k<= i-1;k++){
				y+= a[i][k] * a[k][j];
			}
			a[i][j] -= y;
		}
		imax = j;
		big = fabs( a[j][j] ) / max[j];
		for(i=j+1;i<n;i++){
			y=0;
			for(k=0;k<=j-1;k++){
				y+= a[i][k] * a[k][j];
			}
			a[i][j] -= y;
			if( y = fabs( a[i][j] ) / max[i] > big ){
				big = y;
				imax = i;
			}
		}
		if( imax != j ){
			pivot[j] = imax;
			sign *= -1;
			y = max[j];// exchange max[j] and max[imax]
			max[j] = max[imax];
			max[imax] = y;
			for(k=0;k<n;k++){// exchange line imax and line j
				y = a[j][k];
				a[j][k] = a[ imax ][k];
				a[ imax ][k] = y;
			}
		}
		else{
			pivot[j] = j;
		}
		for(i=j+1;i<n;i++) a[i][j] /= a[j][j];
		//for test only
		/*
		cout<<"~~~~~~~~~ j="<<j<<endl;
		for(i=0;i<n;i++){
			for(k=0;k<n;k++){
				cout<<a[i][k]<<",";
			}
			cout<<endl;
		}
		*/
	}
	delete [] max;
}

/*
 * ludcmp_solve:  solves linear equation Ax=b, the answer x is written into b
 * a is not A, a is LU decomposition of pivotted version of A, and the pivotting is recorded in pivot
 * When U is singular, throws out an error and exit the program.
 */
void ludcmp_solve(int n, double ** a, int * pivot, double * b){

	int j,k;
	double y;

	// do the pivotting on b
	for(j=0;j<n-1;j++){
		k = pivot[j];
		y = b[j];
		b[j] = b[k];
		b[k] = y;
	}
	//cout<<"b after pivoting: "; for(j=0;j<n;j++)cout<<b[j]<<","; cout<<endl;
	// calculate y = Ux, and write it onto b
	for(j=0;j<n;j++){
		y=b[j];
		for(k=0;k<=j-1;k++){
			y -= a[j][k] * b[k];
		}
		b[j] = y;
	}
	//cout<<"y:"; for(j=0;j<n;j++)cout<<b[j]<<","; cout<<endl;
	// calculate x, and write it onto b
	for(j=n-1;j>=0;j--){
		y=b[j];
		for(k=j+1;k<n;k++){
			y -= a[j][k] * b[k];
		}
		if( fabs(a[j][j])<1E-9 ){
			cout<<"error: U is singular in the LU decomposition of A.\n";
			exit(1);
		}
		b[j] = y/a[j][j];
	}
}

/*
 * TridiagLinear:       Solve linear equations Ax=r, where A has a tridiagonal form:
 *                      b0      c0
 *                      a1      b1      c1
 *                      0       a2      b2      c2
 *                      ...     ...     ...     ...
 *                      ...     ...     ...     ...
 *                                      a_{n-1} b_{n-1}
 *
 *                      The solution is stored in u
 */
void TridiagLinear( int n, double *a, double *b, double *c, double *u, double *r ){

        int j;
        double *gamma = new double [n];
        u[0] = r[0]/b[0];
        gamma[0] = 0;
        for(j=1;j<n;j++){
                gamma[j] = c[j-1]/(b[j-1] - a[j-1]*gamma[j-1]);
                u[j] = ( r[j] - a[j]*u[j-1] )/(b[j] - a[j]*gamma[j]);
        }
        for(j=n-2;j>=0;j--){
                u[j] -= gamma[j+1]*u[j+1];
        }
        delete [] gamma;
}

/*
 * jacobi(...): diagonalizes a symmetric real matrix, outputs eigenvalues and eigenkets
 *
 * n is the dimension of the symmetric matrix
 * a[n*n] is the matrix to be diagonalized, the element Aij=a[i*n+j]
 * z gives the eigenvalues of this matrix A,from small to large
 * v (its column arrays) gives the eigenfunctions corresponding to z. 
 * z's and v's initial values do not matter, they would be assigned in this code
 * e.g. bj=v[j][p](p=0,..,n-1) is the jth(j=0,..,n-1) normalized eigenfuction of z[j].
 */
void jacobi(int n, double *a, double *z, double **v)
{
	int l, i, j, p, q;
	double max, temp, s, c, omega, x, y;
	double g, h;
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			v[i][j] = 0;
		}
		v[i][i] = 1.0;         //v is a unit matrix
	}

	for (l = 1; l < max_iteration; l++)
	{
		max = 0;// find the largest matrix element
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				if ( i != j && fabs( a[i*n+j] ) >max )
				{
					max = fabs( a[i*n+j] );
					p = i;
					q = j;
				}
			}
		}

		if ( max < 1e-10 ) 
			break;

		x = -a[p*n+q];
		y = ( a[q*n+q] - a[p*n+p] )/2.0; 
		omega = x / sqrt(x*x+y*y);
		if (y < 0) omega = -omega;
		s = 1.0 + sqrt(1.0 - omega*omega);
		s = omega / sqrt(2.0*s); 
		c = sqrt(1.0 - s*s);
		temp = a[p*n+p];
		a[p*n+p] = temp * c * c + a[q*n+q] *s *s + a[p*n+q] *omega;
		a[q*n+q] = temp * s * s + a[q*n+q] *c *c - a[p*n+q] *omega;
		a[p*n+q] = 0;
		a[q*n+p] = 0;

		for (j = 0; j < n; j++)
		{
			if (j != p && j != q)
			{
				temp = a[p*n+j]; 
				a[p*n+j] = temp * c + a[q*n+j] * s;
				a[q*n+j] = -temp * s + a[q*n+j] * c;
			}
		}
		for (i = 0; i < n; i++)
		{
			if (i != p && i != q) 
			{
				temp = a[i*n+p];
				a[i*n+p] =  temp * c + a[i*n+q] * s;
				a[i*n+q] = - temp * s + a[i*n+q] * c;
			}
		}
		for (i = 0; i < n; i++)
		{
			temp = v[i][p];
			v[i][p] =  temp * c + v[i][q] * s;
			v[i][q] = - temp * s + v[i][q] * c; 
		}
	}//l
	if( l== max_iteration ){
		cout<<" Jacobi: after "<<max_iteration<<" iterations, didn't converge...\n";
		exit(1);
	}
	for (i = 0; i < n; i++)
		z[i] = a[i*n+i];
	for(p=0;p<n-1;p++)
	{
		temp=z[p];
		l=p;
		for(i=p+1;i<n;i++){// bubble sorting
			if(g>z[i]){// find the smallest eigenvalue in z[p], ..., z[n-1]
				temp=z[i];
				l=i;
			}
		}
		if(l!=p)
		{
			temp=z[p];// put the smallest eigenvalue in {z[p],...,z[n-1]} in z[p]
			z[p]=z[l];
			z[l]=temp;
			for(j=0;j<n;j++)// exchange v[][p] with v[][l]
				// this can be optimized, we can store eigenwave in a row, rather than in a column
			{
				temp = v[j][p];
				v[j][p]=v[j][l];
				v[j][l]=temp;
			}
		}
	}                      //eigvalues from small to large
	for(i=0;i<n;i++)
	{
		temp=0;
		for(j=0;j<n;j++)
		{
			temp+=v[j][i]*v[j][i];
		}
		for(j=0;j<n;j++)
		{
			v[j][i] /= sqrt(temp);  //normalization of wavefunctions
		}
	}
}
